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Based on the Suzuki product-formula approach, we construct a family of unconditionally stable 
^ ' algorithms to solve the time-dependent Maxwell equations. We describe a practical implementation 

of these algorithms for one-, two-, and three-dimensional systems with spatially varying permittivity 
04 ' and permeability. The salient features of the algorithms are illustrated by computing selected eigen- 

modes and the full density of states of one-, two-, and three-dimensional models and by simulating 
the propagation of light in slabs of photonic band-gap materials. 

^ ' PACS numbers: 02.60.Cb, 03.50.De, 41.20.Jb, 42.70.Qs 



I. INTRODUCTION 



Maxwell's partial differential equations of electromagnetism describe the evolution of electric and magnetic fields 

in time and space W . They apply to a wide range of different physical situations that are specified by the boundary 

Q , conditions on the electromagnetic (EM) fields. In many cases, numerical methods are required to solve Maxwell's 

■ rr^ ■ equations, either in the frequency or time donrain. For the time domain, a well-known class of algorithms is based on 

j>^' a method proposed by Yee |0| and is called the finite-difference time-domain (FDTD) method. The FDTD method 

(~| has matured during past years and various algorithms implemented for specific purposes are available by now [^,Q. 

CIh These algorithms owe their popularity mainly due to their fiexibility and speed while at the same time they are easy 

'~~', to implement. A limitation of Yee-based FDTD techniques is that the stability of the algorithms is conditional. The 

,__! ■ stability depends on the mesh size used to discretize space and the time step used to perform the time integration [|j . 

J>. , In this paper we present a family of unconditionally stable algorithms that solve the time-dependent Maxwell 

C<1 ' equations (TDME) through the application of orthogonal transformations. That this is possible follows from the 

O^ [ representation of the TDME in matrix form. The exponential of a skew-symmetric matrix plays the role of the time- 

^^ ■ evolution operator of the EM fields. This time-evolution operator is orthogonal. The key to the construction of an 

unconditionally stable algorithm to solve the Maxwell equations is the observation that orthogonal approximations 

to this operator automatically yield unconditionally stable algorithms. The Lie- Trotter-Suzuki product formulae |5[| 

provide the mathematical framework to construct orthogonal approximations to the time-evolution operator of the 

Maxwell equations. However, this framework does not specify how to implement the algorithm. 

Q ' Recently, a spectral-domain split-operator technique has been proposed to solve the TDME Q . The split-operator 

approach is based on one of the many forms of the Lie- Trotter-Suzuki product formulae. The spectral-domain 

^1^' method makes use of Fast Fourier Transforms to compute the matrix exponentials of the displacement operators. 

r^ I The choice made in Ref. [|6| yields an approximation to the time-evolution operator that is no longer orthogonal and 

rS |, hence unconditional stability is not automatically guarantueed M. In contrast, the methodology that we propose 

IJ ' yields efficient, explicit, unconditionally stable schemes that operate on the EM fields defined on the real space grid 

• ^H , only. This renders the algorithms rather ffexible, avoids wrap-around effects H , and naturally allows for the spatial 

/\^ • variations in both the permittivity and the permeability. On the other hand, the implementation described in this 

JH , paper is by no means unique and leaves a lot of room for further improvements. 

■ - - ■ For EM fields in a homogeneous medium, Zheng et al. showed that there is an alternating-direction-implicit (ADI) 



time-stepping algorithm that is unconditionally stable |8||9(]. Conceptually this approach is different from ours. The 
Fourier-mode stability analysis performed by Zheng et al. ||] does not generalize to the case of spatially varying 
permittivity and permeability whereas in our approach, the algorithms are unconditionally stable by construction. 
Our presentation is organized as follows: The basic theoretical concepts are given in Sec. 0. In Sec. HI we explain 



the general philosophy that underlies the Suzuki approach to construct algorithms that are unconditionally stable. 



In Sec. IV , we show in detail how to implement these ideas for the case of the TDME in one, two and three spatial 
dinrensions, using algorithms that are accurate up to second and fourth order in the time step. In Sec. M we explain 



how we analyze the data generated by these algorithms. In Sec. VI we present the results of numerical simulations 



for the physical systems that we selected as examples to test the algorithms. Our conclusions are given in Sec. VIl 
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II. THEORY 

The model system we consider in this paper describes EM fields in a rf-dimensional {d = 1,2,3) medium with 
spatially varying permittivity and/or permeability, surrounded by a perfectly conducting box. In the absence of free 
charges and currents, the EM fields in such a system satisfy Maxwell's equations ||l| 

(1) 

(2) 

(3) 

(4) 

where H = {Hx{r,t),Hy{r,t),Hz(r,t))'^ and E = {Ex{r,t), Ey{r,t),Ez{r,t))'^ denote the magnetic and electric field 
vector respectively. The permeability and the permittivity are given by ii — /i(r) and e — e(r). For simplicity of 
notation, we will omit the spatial dependence on r = {x, y, z)^ unless this leads to ambiguities. On the surface of the 
perfectly conducting box the EM fields satisfy the boundary conditions [y 

nxE = , nH = 0, (5) 

with n denoting the vector normal to a boundary of the surface. The conditions Eq. (ph assure that the normal 
component of the magnetic field and the tangential components of the electric field vanish at the boundary ||I| . Some 
important symmetries of the Maxwell equations (|l|)-(§) can be made explicit by introducing the fields 

X(t) = VMH(i) , Y(t)-ViE(f). (6) 

In terms of the fields X(i) and Y(i), the TDME (Eqs.(|i|) and (|)) read 

a^ U(t) j - 1^ -ijV X -^ )\Y{i))-^{Y{t))- (^) 

Writing m{t) = (X(t),Y(i))^, Eq. (0) becomes 

^^^[t)^n^{t). (8) 

It is easy to show that Ti is skew-symmetric, i.e. 7i^ = — 7i, with respect to the inner product 

(*!*') = / *^-*'dr, (9) 

Jv 

where V denotes the volume of the enclosing box. 
The formal solution of Eq. (||) is given by 

*(i) = e*^*(0), (10) 

where \I'(0) represents the initial state of the EM fields and the operator 

U{t) = e*«, (11) 

determines their time evolution. By construction 

\\^{t)f = {-^{mii)) = I [e^\t)+^^n\t)] dv, (12) 

Jv 

relating the length of '^{t) to the energy density 

w;(i)=£E2(i)+/iH2(i), (13) 

of the EM fields §. From U{tY = U{-t) = U-^{i) = e"*^ it follows that (C/(t)^(0)|C/(i)1'(0)) = (*(i)|*(t)) = 
(^(0)1^(0)). Hence the time-evolution operator U{t) is an orthogonal transformation, rotating the vector ^(t) without 



changing its length |1^||. In physical terms this means that the energy density of the EM fields does not change with 
time, as expected on physical grounds [y. 

The fact that U{t) is an orthogonal transformation is essential for the development of an unconditionally stable 
algorithm to solve the Maxwell equations. In practice, a numerical procedure solves the TDME by making use of an 
approximation U{t) to the true time evolution U{t) (see below). A necessary and sufficient condition for an algorithm 
to be unconditionally stable is that (lO| 

||(7(i)*(0)||<||*(0)||. (14) 

In other words, the length of ^(t) should be bounded, for arbitrary initial condition ^/{t = 0) and for any time 
t [Q. By chosing for ^(0) the eigenvector of U{t) that corresponds to the largest eigenvalue of U{t), it follows from 
Eq. (nj) that the algorithm will be unconditionally stable by construction if and only if the largest eigenvalue of U{t) 
(denoted by ||C/(t)||) is less or equal than one [|lO|. If the approximation U{t) is itself an orthogonal transformation, 
then II C/(t) II — 1 and the numerical scheme will be unconditionally stable. 

Summarizing: Unconditionally stable algorithms to solve Eq. (M) can be constructed by employing orthogonal 
approximations to the time-evolution U{t) = e*'^. For the case at hand, unconditional stability is tantamount to the 
exact conservation of the energy density. 

III. UNCONDITIONALLY STABLE ALGORITHMS 

A numerical procedure that solves the TDME necessarily starts by discretizing the spatial derivatives (see Secjrv|). 
This maps the continuum problem described by Ti. onto a lattice problem defined by a matrix H. Ideally, this mapping 
should not change the basic symmetries of the original problem. The underlying symmetry of the TDME suggests 
to use matrices H that are real and skew-symmetric. Formally the time evolution of the EM fields on the lattice is 
given by 

*(i + T) = C/(T)*(i) ==e^'^*(t). (15) 

The second ingredient of the numerical procedure is to choose an approximation of the time-step operator U{t). A 
standard procedure is to truncate the Taylor series of the matrix exponential |pl],p^ 



U{r) =e^" = f: ^. (16) 

^ — ' n! 

n=0 

Retaining terms up to first order in t yields 

U{t)=I + tH, (17) 

where / denotes the identity operator. As U{t)U{t)'^ — I — (tH)'^ ^ I ior t j^ 0, it is clear that the matrix 
Eq. dl^) is not orthogonal. Making use of the symmetry of H and the positivity of the inner product, we find that 
(*(r)|*(T)) = (c7(r)1'(0)|f7(T)*(0)> ^ (*(0)|*(0)) -h T2(7J*(0)|iJ*(0)) > (*(0)|*(0)), implying that Eq. (0) does 
not hold. Hence, according to the arguments given above, the (Euler) scheme Eq. (O) is unstable, a fact which is of 
course well-known Q . 

The Yee algorithm is based on a leapfrog arrangement pj and formally corresponds to an approximation to the 
matrix exponential U{t) that can be written as 

UYee{T) = I + THi+T^H2, (18) 

where Hi and H2 are matrices, the structure of which depend on the lattice dimensionality. The presence of the 
second order contribution can render the algorithm stable under certain conditions. It seems difficult to determine 
these conditions for arbitrary (skew-symmetric) Hi and (symmetric) H2, i.e. without making use of very specific 
knowledge about the elements of Hi and H2. For EM fields moving in free space, a Fourier-space stability analysis of 
the Yee algorithm yields 

as the condition for stability H . Here c is the light velocity in vacuum and A denotes the spatial mesh size H . 




A systematic approach to construct orthogonal approximations to matrix exponentials, i.e. to construct uncondi- 
tionally stable algorithms, is to make use of the Lie- Trotter-Suzuki formula ||l^,14| 

(20) 

and generalizations thereof ||,|l^. Applied to the case of interest here, the success of this approach relies on the basic 
but rather trivial premise that the matrix H can be written as 

p 
H = J2h„ (21) 

where each of the matrices Hi is real and skew-symmetric. 



The expression Eq. (20) suggests that 



[/i(r) =e^^^..e^«^ (22) 

might be a good approximation to U{t) if r is sufficiently small. Most importantly, if all the H^ are real and skew- 
symmetric, Ui{t) is orthogonal by construction. Therefore, by construction, a numerical scheme based on Eq. (p3) 
will be unconditionally stable. Using the fact that both U{t) and Ui{t) are orthogonal matrices, it can be shown 
that 

||C/(r)-C/i(r)||<y^||[i/„if,]||, (23) 

i<j 

where [Hi.Hj] = HiHj — HjHi. From Eq. ( ^ ) it follows that, in general, the Taylor series of U{t) and Ui{t) are 
identical up to first order in r. We will call Ui\t) a first-order approximation to U{t). 

The product-formula approach provides simple, systematic procedures to improve the accuracy of the approximation 
to U{t) without changing its fundamental symmetries. For example the orthogonal matrix 

U2{t) = [/i(-r/2)^C/i(T/2) ^e^"^/^ ...e^"^/^e^"^/^ ...e-""^/^, (24) 

is a second-order approximation to U{t) |a,^9|. Suzuki's fractal decomposition approach [H gives a general method 



to construct higher-order approximations 
is given by g 



jased on Ui(t) or U^ij). A particularly useful fourth-order approximation 



Ui[T) = C/2(aT)C/2(aT)C/2((l - 4a)T)C/2(aT)t/2(aT), (25) 



where a = ^IjA ~ 4^/'^). The approximations Eqs.(p2|) and (^J), and ( psj) have proven to be very useful in many 
applications J14M23] and, as we show below, turn out to be equally useful for solving the TDME. In practice an efficient 
implementation of the first-order scheme is all that is needed to construct the higher-order algorithms Eqs.(p4[) and 

To summarize: Suzuki's product-formula approach provides the formal machinery to define algorithms that are 
unconditionally stable by construction. The accuracy of these algorithms can be improved systematically, to any 
desired order B. The only assumption made so far is that the real, skew-symmetric matrix H representing the 
TDME can be written as a sum of p real, skew-symmetric matrices Hi. The next step is to choose the HiS such that 
the matrix exponentials exp(r-ffi), ..., exp(r_ffp) can be calculated efficiently. This will turn the formal expressions 
for U2{t) and ^^(t) into efficient algorithms to solve the TDME. 

IV. IMPLEMENTATION 

In this section we present the details of our implementation of unconditionally stable algorithms to solve the 
TDME based on the Suzuki product-formula approach. For pedagogical reasons we start by considering the simplest 
case: A one-dimensional (ID) system. Then we show that the strategy adopted for ID readily extends to higher 
spatial dimensions. The implementation we describe below is by no means unique, leaving a lot of room for further 
improvements. In principle any decomposition Eq. (El]) of H into real skew-symmetric parts will do. Largely guided 



by previous work |17| , |l9| , p4y25[ , we have adopted a decomposition that is efficient, flexible, sufficiently accurate and 
easy to program. 



A. One dimension 

We consider a ID system along the x-direction. Accordingly, Maxwell's equations contain no partial derivatives with 
respect to y or z and e and ^ do not depend on y or z. Under these conditions, the TDME reduce to two independent 
sets of first-order differential equations M . The solutions to these sets are known as the transverse electric (TE) mode 
and the transverse magnetic (TM) mode ||l| . As the equations of the TE- and TM-mode only differ by a sign we can 
restrict our considerations to the TM-mode and obtain the result for the TE-mode by reversing the time. 

From Eq. (Q) it follows that the magnetic field Hy{x,t) = Xy{x^t)/ yj ^{x) and the electric field Ez{x,t) = 
Yz{x, t)/ ^Je{x) of the TM-mode are solutions of 

^X,(..t)^^Af^V (26) 

dt y^ dx y yTJ^Xx) J 

Note that in ID, the divergence of Hy{x,t) and Ez{x,t) is zero. Hence Eqs.(0) and (0) are automatically satisfied. 
Using the second-order central-difference approximation to the first derivative with respect to x, we obtain 

|x.(. t) ^ -^ f ^:^(i±M _ M_AA) , (28) 

dt 5^ V ^fWJ+l n/T^TT / 

where the integer i labels the grid points and 5 denotes the distance between two next-nearest neighbor lattice points 
(hence the absence of a factor two in the nominator). For notational simplicity we will, from now on, specify the 
spatial coordinates through the lattice index i, e.g. Xy{i^ t) stands for Xy{x = (i + l)(5/2, t). 

/ Hy E^ Hy E^ Hy E^ Hy IX 

/ •— — • •— — •- - - - - •— • •^X 



/] 1 2 3 4 n-2 n-1 n / 

8 

FIG. 1. Positions of the two TM-mode field components on the one-dimensional grid. 

Following Yee JH| it is convenient to assign Xy{i^ t) and Yz{j, t) to the odd, respectively, even numbered lattice site, 
as shown in Fig. [Jfor a grid of n points. The equations (^8|) and ( p9| ) can now be combined into one equation of the 
form Eq. (|8|) by introducing the n-dimensional vector 

*(«, t) = I J^/,*'/) = V^f y(*' ^)' \ «dd 

\Yz{i,t) = ^iEy{i,t), I even ^ ' 

The vector ^(t) describes both the magnetic and the electric field on the lattice points i — 1, . . . ,n. As usual the 
i-th element of ^(i) is given by the inner product ^(i,t) = ef ■ ^(i) where e^ denotes the z-th unit vector in the 
n-dimensional vector space. Using this notation, it is easy to show that Eqs.(p8[) and ( |29| ) reduce to 



|*(i) = H*(i), (31) 

where the matrix H is given by 



^ = Y1 t'^^+i.* i^i •^^+1 " ^^+l^f) + /3"+i,»+2 (ei+ief+2 - e,+2ef+i)] , (32) 



E 



with (3ij — l/{S^/eiJIj) and the prime indicates that the sum is over odd integers only. In complete analogy to Eq. ( [lO| ) 
the time evolution of ^(i) is formally given by *(i) = U{t)^{0) with U{t) = cxp{tH). 



The notation introduced above will prove most useful for the case of 2D and 3D for which it is rather cumbersome 
to write down matrix representations. For the ID case it is not difficult and in fact very instructive to write down 
the matrix H explicitly. Indeed, we have 



H = 



( /?2,1 

-/32,i 



V 



-/?n-l,n-2 



\ 



I3n~l,7i 



(33) 



and we immediately see that H is ske"w-sym.metric by construction. Furthermore, for n odd we have 

— *(!,<) =/32,i*(2,i) and — *(n, t) = -/3„_i.„*(n - l,i) , (34) 

ot at 

such that the electric field vanishes at the boundaries (^2(0, t) =Yz{n + 1, t) ~ 0, see also Fig. |l|), as required by the 
boundary conditions Eq. dq). For this reason we only consider the case of n odd in the sequel. 

According to the general procedure outlined in Sec. Ill, the final step in the construction of an unconditionally stable 
algorithm is to decompose H. Guided by previous work on Schrodinger and diffusion problems [|l7| , [l9| , p4|j25|] , we split 
H into two parts, i.e. H = Hi + H2, where 



n 

H2 = ^ A+i,j+2 (e^+ie 



■'1+2 



i+2^i+lj I 



(35) 
(36) 



i.e. we divide the lattice into odd and even numbered cells. In matrix notation we have 



/ /32,i 

-/32,i 





Hi 





/34,3 

-/34,3 





Pn- 



l,n-2 







0/ 



(37) 



and 







Ho 





/32,3 

-/32.3 





-Pn- 



V 








/?„_3.„_2 




-3,n-2 








/?„_!,„ 




~/3„-i,„ 



(38) 



Clearly both Hi and H2 are skew-symmetric block-diagonal matrices, containing one 1x1 matrix and [n— l)/2 real, 
2x2 skew-symmetric matrices. 

According to the general theory given above, the first-order algorithm defined by 



Ui{t) 



tHi TH2 



(39) 



is all that is needed to construct unconditionally stable second and higher-order algorithms. As the matrix exponential 
of a block-diagonal matrix is equal to the block-diagonal matrix of the matrix exponentials of the individual blocks, 
the numerical calculation of e^^^ (or e^^^) reduces to the calculation of (n— 1)/2 matrix exponentials of 2 x 2 matrices. 
The matrix exponential of a typical 2x2 matrix appearing in e^^^ or e^^^ is given by 



exp 



1 
-1 



*(*,t) 
*(i,t) 



cos a 
- sin a 



smck 
cos a 



*(j,t) 



(40) 



B. Two dimensions 



Assuming translational invariance with respect to the z direction, the system effectively becomes 2D and the TDME 
separate into two sets of equations ||l[. For conciseness, in this section we only discuss the set of equations for the 
TM-modes. The TE-modes can be treated in exactly the same manner. 

The relevant EM fields for the TM-modes in 2D are \E'(i) = (Xa;(x, ?/,f),Xj,(x, y, t), Y2(a;,y, t))"^, in terms of which 
TDME Eq.(0) reads 



d_ 

di 



*(i) = W*(t) 



j__a_j_ 1 d 1 




(41) 



We discretize continuum space by simply re-using the one-dimensional lattice introduced above, as exemplified in 
Fig. S for the case of the TM-modes. This construction automatically takes care of the boundary conditions if n^ and 
Uy are odd and yields a real skew-symmetric matrix H. In analogy with the ID case the elements of ^(i) are defined 
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FIG. 2. Positions of tlie three TM-mode EM-field components on the two-dimensional grid for Ux = 9 and Uy — 5. 



Yz{i,j,t) — y/eTJEz{i,j,t), i even and j even 
^ihJit) = { Xy{i,j,t) = ^/JI~JHy{i,j,t), i odd and j even 
Xx{i,j,t) = ^/JI~JHx{i,j,t), i and j odd 



(42) 



Discretization of the differential operators that appear in Eq. (B^) yield expressions that have the same structure as 
Eq. (|33|), with extra subscripts to account for the second spatial dimension. It follows that on the lattice. 



d_ 
dt 



*(i) = H^{t) - Y.'T.' [H^^Hhj) + H^''Hhj)] *(i), 



1=1 J=l 



where 



j^(x)/j -N ^ I ^i,j+l^i+l,j + l ^i+l,] + l^i.] + l ^t+l,] + l^i+2,] + l ^■i+2,j + l^i+l,] + l 



H^y\i,j)^ 



T T 

S ^/el+TJ+TjH+TJ 



T T 



(43) 

(44) 
(45) 



and the superscripts (a;) and (y) refer to the derivative with respect to x and y respectively. Note that we use the 
pair (i,j) to label the UxUy unit vectors etj. In complete analogy with the ID case we split Eqs.([44|) and ( [45| ) into 
two parts and obtain for the first-order approximation to U{t) 



Ui{t) 



rfff"' „rH<"> „rH<''' „rH<''' 



(46) 



where for instance, in formal analogy to Eq. (36), we have 






' ^i+lj + l^i+2,] + l 



^i+2.] + l^i+lj + l 



Sy^£i+2j+TjM+2~j + T 



(47) 



It is not difficult to convince oneself that approximation Ui{t) and hence also U2{t) and U/^{t) do not commute with 
the (lattice version of) the divergence. Therefore the divergence of the EM fields in 2D is not conserved. Although 
the initial state (at i = 0) of the EM fields satifies Eqs.(p|) and (0), time-integration of the TDME using Ukij) yields 
a solution that will not satisfy Eqs.(^) and (Q). However, for algorithm Ukir) the deviations from zero vanish as r*"' 
so that in practice these errors are under control and can be made sufficiently small for practical purposes. 



C. Three dimensions 



In terms of *(t) = (X(t), Y(i))^ for a 3D system, Eq. (0) reads 



-vl'(i) = n^{t) 
where h is given by 



h 
-K^ 



*W, 



/ 







h = 



J__9_J_ ]_d_J_ \ 



J_^J_ 13 1 





^ dz ^ 

I J_9_J_ 1 d 1 

\ yTi- 9y ^ yTJI dx ^ 



(48) 



(49) 



We discretize the spatial coordinates by adopting the standard Yee grid |^ . We show a unit cell of this grid in Fig. ||. 



ZA 



y 



FIG. 3. Unit cell of the Yee grid. 



In analogy to the systems in ID and 2D, we assign the EM fields to the lattice points such that the boundary conditions 
Eq. (p) are automatically satisfied. The elements of the vector '^{i,j,k,t) are given by 



^{i,j,k,t) 



Xx{i,j,k,t) = ^/JI~JJ^Hx{i,j,k,t), i even, j odd, k odd 

Xy{i,j,k,t) = ^fiij^kHy{i,j, k, t), i odd, j even, k odd 

Xz{i,j,k,t) = ^JjMjj:Hz{i,j,k,t), i odd, j odd, k even 

Yx{i, j, fc, t) = y^sTJj^Ex{i, j, fc, t), i odd, j even, k even 

Yy{i,j,k,t) = y^EiJj^Eyii, j,k,t), i even, j odd, k even 

Yz{i,j,k,t) = y/eiJ^Ez{i,j,k,t), i even, j even, k odd 



(50) 



for the origin of the coordinate system {i,j, k) ~ (0,0,0) at the center of the unit cell shown in Fig. |[ The number 



of lattice points in the x, y, and z direction will be denoted by rix 
are assumed to be odd. 



'y. 



and Hz respectively. As before these numbers 



Discretization of the differential operators that appear in Eq. ( [49| ) yields Eq. (0) in the form 

n^ ny rij. 

^vI/(i)^i/vI,(i)^^'^'^'[i7(-)(^,J,A;)+i?(^)(^,J■fc) + H(^)(^,J■fc)]vI/(i), (51) 

i=i j=i fc=i 

where the superscripts (x), (?/) and (z) refer to the derivative with respect to x, y and z respectively, 

T T T T 

'^•\/^i+l,i + l,feMij + l,fc '^•\/£i+l,j,fc+l/^ij\fc+l 

/ji jn 2^ 2^ 

^i+l,j + l,k^i+2,j + l,k ^ ^i+2,j + l,k^i+l,j + l,k ^i+l,j,k+l^i+2,j,k+l ^ ^i+2,j,k+l^i+l,j.,k+l 



S y/£i+Tj+Tjrpi+2j+T^ S y/Si+rjjt+lTk+2JM^ 



and the expressions for H^y^ijjjk) and H^^\i,j,k) follow from Eq. ( [5^ ) by symmetry. Note that we use the triple 
(i, j, k) to label the n^nynz unit vectors Gij^k- 

In complete analogy with the ID and 2D case we split Eq. (p2|) (as well as H^y\i^j,k) and H'-^\i,j,k)) in two 
parts and obtain for the first-order approximation to U{t) 

U,{r) = e-<'e^^^^'e^<'e^<'e^<'e-<', (53) 

where for instance 

j^{z) ^ sr^' sr^' sr^' I ^i,]+i,k^i,]+i,k+i ^ ^i,j+i,k+i^i,j+i,k _ ^i+i.j,k^i+i,],k+i " ^i+i,],k+i^i+i,j,k \ ,^^-. 

"'i j^i "^1 \ S^£u]+i,k+i^i,j+i,k 6^/eI^T~jI+TjM+T~j^ J 

Note that each contribution to Eq. (|54| ) acts on a different pair of elements of ^(t). Hence each of the matrix 
exponentials in Eq. (M) acts on one quarter of all the lattice points. Performing the time-step operation Eq. ( [5^ ) 
involves only two sweeps over all lattice points. 

By construction the algorithm defined by Eq. (BSl) is unconditionally stable and so are the higher-order algorithms 
defined by U2{t) and ^^(t). Each contribution to e.g. Eq. ( |54| ) is of the form Eq. ( |35|) and hence its matrix exponential 
can be calculated in exactly the same manner as in the ID case (see Eq. (M)). The divergence of the EM fields in 
3D is, for the same reason as in the 2D case, not conserved but decreases as t^ . 

D. Implementation: Summary 

The notation required to write down the algorithms in mathematical form might give the impression that these 
algorithms are difficult to program. Actually that is not the case, on the contrary. Recall that the first-order algorithm 



C/i(r) is all we need to program: As explained in Sec. Ill, more accurate schemes can be implemented without extra 



programming. Let us consider the algorithm for the case of ID. For 2D and 3D we simply repeat the steps described 
below two, respectively, three times. We have 

= < n "^^P [/^»+i>' (^» ^^+1 ^ ^^+l^^)] M n '^^P [/5«+i^«+2 {^i+i^J+2 - e»+2e^i)] \ ■ (55) 

where we used the block-diagonal structure of Hi and H2 (see Eqs.(|37|) and (pq)) to obtain an exact expression for 
Ui{t) in terms of an ordered product of matrix exponentials. Each of these matrix exponentials only operates on a 
pair of elements of ^(i) and leaves other elements intact. The indices of each of these pairs are given by the subscripts 
of e and e^. From Eq.(|5g) it is then clear what the program should do: Make loops over i with stride 2. For each 
i pick a pair of elements from '^{t) according to the subscripts of e and e^, compute (or fetch from memory) the 
elements of the plane rotation (see Eq. (^)), perform the plane rotation, i.e. multiply the 2x2 matrices and the 
vectors of length two, and overwrite the same two elements. 

It also follows immediately that performing a time step with algorithms based on Eq. (^5|) takes 0{K) plane 
rotations where K is the total number of elements of the vector ^(i) (which is less or equal to the number of grid 
points). This renders the algorithm efficient: The number of operations to complete one time step scales linearly with 
the number of grid points. Also note that there is a high degree of intrinsic parallelism in this class of algorithms. In 
principle, the {n — l)/2 matrix-vector multiplications that implement e'^^^ or e'^^'^ can be done in parallel. 



In the absence of external currents (see below), updating the EM field values of a 3D system using the the Yee 
algorithm requires 6 arithmetic operations (see Eq.(33) in Ref. 0) whereas the second-order algorithm U2{t) requires 
33 arithmetic operations. For ID and 2D problems the ratio is 9/4 and 21/6 respectively. Thus, in terms of CPU 
time, the price paid for the unconditional stability of the algorithms is not that much and for some applications (see 
below) may well be worth paying. 

An important aspect that we have not yet discussed is the effect of the discretization of space on the accuracy 
of the numerical results. Both conditional Yee- type algorithms and unconditionally stable algorithms Ui{t), U2{t) 
and U4{t) suffer from numerical dispersion (see Ref. [^, chapter 4 for an in-depth discussion). Simple methods to 
reduce numerical dispersion are taking a finer mesh or employing more accurate finite-difference approximations for 
the spatial derivates Q. The former obviously can be used here too (for the simulations discussed below we used a 
mesh size that yields sufficiently accurate results for the present purposes). There are no fundamental nor practical 
problems to incorporate the latter method in the Suzuki-product-formula approach |T^,n9] . However as the emphasis 
of the present paper is on the construction of unconditionally stable algorithms we relegate a presentation of these 
technical but for applications important extensions to future publications. 

V. DATA ANALYSIS 

Time-domain algorithms obviously yield the time development of the EM fields. The scattering (transmission) 
of the EM fields from (through) objects is one of the main applications of this technique @]. One approach is to 
prepare an initial state ^(0) of the EM fields, propagate the fields in time for a number of time-steps and analyse 
the scattered and/or transmitted fields. Another, more realistic, approach is to use a current source J(i) = J(r,i). 
Instead of Eq. (ph we have 

^^^i,{t) = n^it)-j{t), (56) 

where 'i'{t) = (X(^), Y(i))-^ and J'(t) = {0,3{t))^ represents the source term. The formal solution of Eq. (pfl) is given 

by 

*(i) = e*"^(0) - / e(*-")^^(u) du, (57) 

showing that we can simply re-use one of the unconditionally stable algorithms to compute the second term in Eq. (|57|). 
In practice, for a time-step r, we update ^(i) according to 

rt+T 

^{t + t) = e^^*(i) - / eS'+^-'^^'^Jiu) du. (58) 



A standard quadrature formula can be used to compute the integral over u p4| . When a current source is present we 
take as the initial condition \l/(0) — 0. 

Time-domain algorithms can also be used to compute the eigenvalues of H , the discretized form of Ti. In general H 
is a (very) large matrix, usually too large to be stored in memory. If only a few, well-separated eigenvalues of H are 
required sparse- matrix techniques can be used to compute these eigenvalues ll2,M. However, if one is interested in 
global features of the distribution of eigenvalues, i.e. if we want to determine all eigenvalues, time-domain algorithms 
offer several advantages. In fact they are at the heart of so-called "fast" algorithms to compute the density of states 
(DOS) and other related quantities |2^-^^. The basic idea of this approach was laid out by Alben et al. |£8[ who 
used it to compute the DOS of models for one electron moving in a disordered alloy. 

Denoting the (unknown) eigenvalues and (unknown) eigenvectors of H by lEj and 4)j respectively, we have 

^^ ^ - (vI'(O)lM'(O)) (*(0)|vI/(0)) ^ (*(0)|vl'(0)) ' ^ ^ 

where K {K ~ n for ID, K — SuxTiy/A for 2D, K — 3nxnynz/4: for 3D) is the dimension of the vector space on 
which H acts. From Eq. ( p9| ) it follows immediately that the Fourier transform of f{t) contains the information on 
all eigenvalues for which |(^(O)|0i)| > 0. Using independent random numbers to initialize the elements of ^(0), it 
can be shown that the density of states V^u) is given by pq] 

/-(-oo 

e-"^'f{t)dt, (60) 

-OO 

10 



where a is an irrelevant constant factor and f{t) is the average of f{t) over different realizations of the random initial 
state. It is often expedient to consider, in addition to I'(w), the integrated density of states 

V{u)du. (61) 

-oo 

The statistical error on f{t) vanishes as l/wSK where S is the number of statistically independent samples of ^(0) 
p6| . The fact that the statistical error decreases with the number of lattice points K/2 gives a tremendous boost to 
the efficiency of the method. 

The information on the eigenvalues of H, obtained through the use of a time-domain method is intimately related 
to the unconditional stability of the latter. As f{t) is band-hmited, with frequencies Ej in the interval [— 1|7?||, ||^||] 
(where ||7f || denotes the largest eigenvalue of H in absolute value), it follows from Nyquist's sampling theorem that it 
is sufficient to sample f{t) at regular intervals Ai = 7r/||iJ||. If A^ denotes the number of data points used to sample 
f(t), frequencies Ej that differ less than Ai? — ir/NAt are indistinguishable (although they will all contribute to 
the DOS). Extending the length of the time integration by a factor of two increases the resolution in frequency by 
a factor of two. This is a rather efficient and flexible procedure to trade accuracy for computational resources. One 
may object that by integrating the TDME over longer and longer times (larger and larger N) the error on '^{t) will 
increase, possibly leading to no gain in accuracy at all. However, it has been shown rigorously pq | that the error on 
the eigenvalues of H vanishes as t^ /N if one uses unconditionally stable algorithms based on the second-order Suzuki 



product- formula. The proof given in Ref. |26| applies to the fourth-order algorithm Ua{t) as well, the exponent of r 
being four instead of two. 

In some cases the underlying differential equations and boundary conditions only specify the solution up to a non- 
zero constant. In the calculation of the DOS, the presence of such a constant contribution shows up as a peak at 
zero frequency. In principle this peak can be removed by modifying the random initial state but as the origin of this 
irrelevant artifact is understood, there is little reason of doing this. 

Summarizing: Solving the TDME by the fc'th-order Suzuki product-formula algorithm Uk{T) guarantees that the 
accuracy with which the eigenvalues of H can be determined vanishes as t^ /N , where N is the number of points used 
to sample f{t) (see Eq. (||)). 

VI. SIMULATION RESULTS 

In this section we present simulation results for several physical systems which we selected as examples to test our 
algorithms. For numerical purposes it is expedient to use dimensionless quantities. We will denote the unit of length 
by A and take the velocity of light in vacuum, c, as the unit of velocity. Then time and frequency are measured in 
units of A/c and c/A, respectively. The permittivity e and permeability ^ are measured in units of their corresponding 
values in vacuum. 

A. One dimension 

Let us first consider an empty, one-dimensional cavity with constant permeability /i = 1 and constant permittivity 
e = 1. The eigenfrequencies for a system of length L are given by 11] 

Trfc 
Wfc = ^ , (62) 

where fc = 0, 1, 2, . . . labels the different EM modes. In Fig. ^ we show the density of states \'D{uj) \ obtained according 
to the procedure described in Sec. 0, using the second-order algorithm C/2(t) (solid line) and the standard Yee 
algorithm |g,^ (dashed line). In both calculations the lattice spacing 5 = 0.1 and the time step r = 0.01. Each 
curve shown in Fig. Q is the average of 5 = 10 statistically independent runs, taking A^ = 16384 samples of f{t) (see 
Eq. ([59|)) at time intervals Ai = 0.1. The peaks in Fig. ^ correspond to the exact frequencies (see Eq. (p^) of the 
ID cavity. The background signal produced by the Yee algorithm [g,|| is at least eight orders of magnitude larger 
than that generated by U2{t). The time-step operator of the Yee algorithm is not an orthogonal matrix and hence its 
eigenvalues do not necessarily lie on the unit circle. This is related to the fact that the Yee algorithm is conditionally 
stable PI and leads to fluctuations in the energy density w(i), as illustrated by the dashed line in Fig. |l and to 
negative values of the Fourier transform of f{t) (which is the reason why Fig. ^ shows |I'(cj)| instead of V(uj)). 

In contrast, the unconditional stability of algorithms based on the Suzuki product-formula implies that 'w{t) is 
constant in time. The solid line in Fig. ^ shows that this is indeed the case. 

11 





10-^ 


^ 


lo" 


9. 


10-* 



10 



10 



10 



3 



to 



FIG. 
Yee alf 



4. Density of states of a one-dimensional cavity of length L = 10. Solid line: U2{t) algorithm; dashed line: standard 
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As a second example we consider a one-dimensional stack of dielectric material, schematically shown in Fig. ^. The 
material indices of refraction, denoted by n\ and ?t.2, give rise to a spatially varying permittivity 



, ^ J Til, a X mod(a + b) < a 
^ ' I r7,2 , if a; mod(a + b) > a 



(63) 



In particular, we consider a structure that is known as the quarter-wave stack and is characterized by the relation 

nia = n2&, (64) 

such that the length of the optical path in the two layers is the same. The density of states exhibits a gap centered 
around the midgap frequency p|,p4[| 



UJQ 



2nia 



(65) 



In Fig. we show the density of states V^ui) as obtained by U2(t). Also for these calculations the lattice spacing 
(5 = 0.1 and the time step r — 0.01. Each curve in Fig. ^ is the average of S* = 100 statistically independent runs, 
taking N = 16384 samples of f(t) (see Eq. (|59|)) at time intervals Ai = 0.1. Note that the system length L and the 
(odd) number of lattice points have to be chosen judisciously, otherwise the spectrum will exhibit artifacts (impurity 
states due to one extra grid point). 
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7. Density of states, as obtained by the U2{r) algorithm, of a quarter- wave stack of length L — 24.9 and with parameters 
712 = 4, a = 0.8, b — 0.2 (see FigJd) as a function of the rescaled frequency ti = uj/luo- 



In Fig. g we show the density of states I?(u;) as obtained by the Yee algorithm. Note that the spectral weight can 
take negative values, an unphysical feature that is a manifestation of the fact that the energy of the EM field is not 
conserved. In Fig. |^ we present the integrated density of states iV(w) for both the U2{t) and Yee algorithm. The 
result of the U2 (t) algorithm is in excellent agreement with the analytical calculation pM . 



B. Two dimensions 



A photonic bandgap (PBG) material prohibits the propagation of EM fields in a range of frequencies that is 
characteristic for its structure |33|. A PBG is called absolute if it exists for any wave vector of the EM fields. The 
most common method used to compute a PBG employs a plane- wave expansion to solve the time independent Maxwell 
equations (see e.g. [3^]). This kind of PBG calculation requires a Fourier transform of the unit cell of the dielectric 
structure which is for simplicity considered to be periodic. 

With our time-domain algorithm the existence of a PBG can be demonstrated with relative ease. It suffices to 
compute the spectrum of such a dielectric structure with a random initial state. If the spectrum is gapless there is 
no need to make additional runs. If there is a signature of a gap, it can be confirmed and refined by making more 
runs. As an example we consider a system consisting of a dielectric material pierced by air-filled cylinders [p^ . The 
geometry is taken to be a square parallelepiped of size L = 45.1 that is infinitely extended in the z-direction and 
hence is effectively two-dimensional. In Fig. O we present the results for the PBGs which we obtained for both the 
transverse magnetic (TM) and transverse electric (TE) modes as a function of the filling fraction. The data have 
been generated by means of the algorithm C/4(t) with a mesh size 6 = 0.1 and a time step r = 0.1. To compute the 
DOS we used N = 32768 samples of f{t) at time intervals At = 0.1. Only a single random initial state for the EM 
fields has been used. Replacing the free-end boundary conditions Eq. (|5|) by periodic boundary conditions (results 
not shown) only leads to minor changes in the locations of PBGs. The results shown in Fig. nfl are in good agreement 
with those presented in Ref. [pTj . 

In Fig. O we study the propagation of time-dependent EM fields through the above described PBG material 
consisting of twelve unit cells. The PBG material is placed in a cavity which contains a point source (located to the 
left of the PBG material) that emits radiation with frequency w. The TDME were solved by the U2{t) algorithm 
with 5 = 0.1 and r = 0.01 in the presence of a current source according to Eq. (|5^). The snapshots show the absolute 
intensity E^ of the TM-mode at i = 102.4. The computed DOS of the PBG material is given in Fig. ^ We used 
the Ui{T) algorithm with 5 = 0.1, r = 0.1, and took N = 32768 samples of f{t) at time intervals Ai — 0.1 in 
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FIG. 8. Density of states, as obtained by the standard Yee algorithm, of a quarter-wave stack of length L — 24.9 and with 
parameters ni — 1, n2 = 4, a = 0.8, 6 — 0.2 (see Fig.H) as a function of the rescaled frequency Cj — lo/ujo- Note the difference 
in the vertical scale between FigsJ?! and pi 
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FIG. 9. Integrated density of states as function of the rescaled frequency di — uj/loo, for the same system as in Figs.M and 
Solid line: U2{r) algorithm; dashed line: standard Yee algorithm; dashed-dotted line: analytical result. 



this computation. The presence or absence of gaps in the DOS leads to quahtative changes in the transmitted (and 
reflected) intensities. Since a gap is present in the DOS at uj — 1.89, radiation with this frequency does not easily 
propagate through the (thin slice of) PBG material. On the other hand, the DOS has no gaps aX oj — 1.50 and 
uj — 2.50, so that propagation of EM fields through the PBG material should be possible, as is indeed confirmed by 
Fig. |1|. 

C. Three dimensions 

We first compute the modes of a simple cubic cavity of size L. The eigenfrequencies are given by |l|| 

iJkhn = TTL-Wk^+l^+m^ , (66) 

where fc, Z, and m are non-negative integers. Eigenmodes corresponding to (fc,0,0), (0,^,0) and (0,0, m) are incom- 
patible with the boundary conditions Eq. (0). In Tab. | we present the frequencies of the five lowest eigenmodes of a 
cubic cavity for L = 5. The agreement with the theoretical values is satisfactory. Note that as the frequency increases, 
the deviation from the exact result increases. This is due to the (second-order) finite-difference representation of the 
spatial derivatives on the grid and is a manifestation of the numerical dispersion mentioned earlier. 
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(fc, I, m) 



l-Okl-n 



Simulation 



1,1,0 
1,1,1 
2,1,0 
2,1,1 
2,2,0 



0.889 
1.088 
1.405 
1.539 
1.777 



0.889 
1.089 
1.404 
1.534 
1.771 



TABLE I. Frequencies of the eigenmodes of a cubic cavity of size L = 5. Simulation; values determined from DOS data 
generated using U2{t) with parameter values 5 — 0.2, r — 0.01, A'' = 4096 and At = 0.1. 

As a second 3D example we consider the emission of the EM radiation from a point source located inside dielectric 
material containing spherical voids. A projection of the material onto the x — y (or y — zorx — z) plane is shown in the 
top panels of Fig.^ A point source is placed inside the center void to mimic an atom or molecule that emits a photon. 
The remaining panels in Fig.O show snapshots of the light intensity after an elapsed time t = 384, for different sizes 
of the voids and different values of the permittivity. If the latter is larger than 5, the images no longer depend on the 
value of the permittivity (results not shown). For e — 1.5 (panels (b) and (B)) the EM field easily propagates through 
the sample and leaves the sample from all sides. This is not the case for e — 5: Radiation can leave the sample only at 
those locations where there is very little or no dielectric material left. In other words, light emerges from the sample 
in well-defined directions only. Clearly, much more work is necessary to study the propagation of EM radiation in 
this system as a function of the material parameters, the system size and the frequency of the emitted photons. 

VII. CONCLUSION 

We have introduced a new family of algorithms to solve the time-dependent Maxwell equations. Salient features of 
these algorithms are: 

• rigorously provable unconditional stability for one-, two- and three-dimensional systems with spatially varying 
permittivity and permeability, 

• the use of a real-space (Yee-like) grid, 

• the order of accuracy in the time step can be systematically increased without affecting the unconditional 
stability (in this paper we limited ourselves to second and fourth-order schemes) 

• the exact conservation of the energy density of the electromagnetic field 

• easy to implement in practice 
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FIG. 11. Snapshot of the intensity E^ at t = 102.4. Dimensions of the system: 30 x 12.1; point source located at (6,6) (see 
Fig.0), emitting radiation with frequency uj. 

We have presented results for the density of states of simple cavities and photonic bandgap materials. We demon- 
strated that mathematical properties of the algorithms are such that they can be used to compute the density of 
states with very good accuracy (limited by the accuracy of the spatial discretization used). We gave some illustrative 
examples of scattering of waves by photonic bandgap systems. These examples also served to show that our algorithms 
reproduce known results. The first feature opens up possibilities for apphcations to left-handed materials [^ 39|. We 
intend to report on this subject in the near future. 

Although we believe there is little room to improve upon the time-integration scheme itself (except for using higher- 
order product-formulae) , for some applications it will be necessary to use a better spatial discretization than the most 
simple one employed in this paper. There is no fundamental problem to extend our approach in this direction and we 
will report on this issue and its impact on the numerical dispersion in a future publication. 

The rigorous unconditional stability of the algorithms that we proposed in this paper is a direct consequence of 
adopting a Suzuki-product-formula approach that preserves the fundamental symmetries of the physical system. In 
view of the generic character of this methodology, the approach pursued in the present paper should be useful for 
constructing unconditionally stable algorithms that solve the equations for e.g. sound, seismic and elastic waves as 
well. 
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FIG. 13. Intensity of the EM-fields in and outside a cubic sample of dielectric material containing 3x3x3 spherical voids, 
placed in an empty cavity. The size of the cavity is 12.1 x 12.1 x 12.1, the size of the sample is 9.1 x 9.1 x 9.1, the source is 
located at (5.6, 5.5, 6), the slices shown are at z = 6, the fourth-order algorithm U4,{t) was used with a time step r = 0.075 and 
mesh size 5 = 0.1. (a): Intensity at i = 0.3. The radius of the empty spheres is 1.4. (b): Intensity at f = 384. Same system as 
in (a). The permittivity of the dielectric material e = 1.5. (c): Intensity at i = 384. Same system as in (b). The permittivity 
of the dielectric material e = 5. (A): Intensity at i = 0.3. The radius of the empty spheres is 1.5. (B): Intensity at i = 384. 
Same system as in (A). The permittivity of the dielectric material e — 1.5. (C): Intensity at i = 384. Same system as in (B). 
The permittivity of the dielectric material e = 5. 
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